subroutine rhs(u,v,u_rhs,v_rhs)
implicit none
real*8::u(:),v(:),u_rhs(:),v_rhs(:)
integer::i
real*8,parameter::c=1.,dx=1.

do i= 3, size(u)-2
    
   u_rhs(i)=v(i)
   v_rhs(i)=c**2*(-1./12.*(u(i+2)+u(i-2))+4./3.*(u(i+1)+u(i-1))-5./2.*u(i))/(dx**2)
end do

end subroutine

subroutine extrapolation_bcs(u)
implicit none
real*8::u(:)
integer::n
n=size(u)
u(2)=3.*u(3)-3.*u(4)+u(5)
u(1)=3.*u(2)-3.*u(3)+u(4)
u(n-1)=3.*u(n-2)-3.*u(n-3)+u(n-4)
u(n)=3.*u(n-1)-3.*u(n-2)+u(n-3)

end subroutine

subroutine bcs(u)
implicit none
real*8::u(:)
integer::n
n=size(u)
u(2)=0
u(n-1)=0
u(1)=u(n-2)
u(n)=u(3)
end subroutine


program main
implicit none
integer,parameter::n=83,m=1500
real*8:: u(n,m)=0,v(n,m)=0,u_k1(n)=0,v_k1(n)=0,u_k2(n)=0,v_k2(n)=0,u_k3(n)=0,v_k3(n)=0,u_k4(n)=0,v_k4(n)=0,u_tmp(n)=0,v_tmp(n)=0
integer::i,j
real,parameter::dt=0.25

interface
  subroutine rhs(u,v,u_rhs,v_rhs)
  implicit none
  real*8::u(:),v(:),u_rhs(:),v_rhs(:)
  integer::i 
  end subroutine

  subroutine extrapolation_bcs(u)
  implicit none
  real*8::u(:)
  integer::n
  end subroutine

  subroutine bcs(u)
  implicit none
  real*8::u(:)
  integer::n
  end subroutine


end interface

!initial data
do i=1,n

  u(i,1)=sin(3.1415926*(i-2)/(n-3)*2)

end do

!runge kutta
do i=1,m-1

  call rhs(u(:,i),v(:,i),u_k1,v_k1)
  do j=3,n-2
    u_tmp(j)=u(j,i)+0.5*dt*u_k1(j)
    v_tmp(j)=v(j,i)+0.5*dt*v_k1(j)
  end do
  call  bcs(u_tmp)
  call  extrapolation_bcs(v_tmp) 
  
  call rhs(u_tmp(:),v_tmp(:),u_k2,v_k2)
  do j=3,n-2
    u_tmp(j)=u_tmp(j)+0.5*dt*u_k2(j)
    v_tmp(j)=v_tmp(j)+0.5*dt*v_k2(j)
  end do
  call bcs(u_tmp)
  call extrapolation_bcs(v_tmp)

  call rhs(u_tmp(:),v_tmp(:),u_k3,v_k3)
  do j=3,n-2
    u_tmp(j)=u_tmp(j)+dt*u_k3(j)
    v_tmp(j)=v_tmp(j)+dt*v_k3(j)
  end do
  call bcs(u_tmp)
  call extrapolation_bcs(v_tmp)

  call rhs(u_tmp(:),v_tmp(:),u_k4,v_k4)

  do j=3,n-2
     u(j,i+1)=u(j,i)+(1./6.)*dt*(u_k1(j)+2*u_k2(j)+2*u_k3(j)+u_k4(j))
     v(j,i+1)=v(j,i)+(1./6.)*dt*(v_k1(j)+2*v_k2(j)+2*v_k3(j)+v_k4(j))
  end do
  call bcs(u(:,i+1))
  call extrapolation_bcs(v(:,i+1))

end do

open(file='a.txt',unit=10,action="write")
do i=1,m
 write(10,"(83f9.5)")u(:,i)
end do

end program
